/*
 * =====================================================================================
 *
 *       Filename:  C_T_P2_2D.h
 *
 *    Description:  define the base functions in 2D case
 *        Version:  1.0
 *        Created:  2021/04/03 03时02分41秒
 *       Revision:  none
 *       Compiler:  gcc
 *         Author:  hhxie@lsec.cc.ac.cn
 *        Company:  
* =====================================================================================
 */
#ifndef __CTP22D__
#define __CTP22D__

#include <stdbool.h>
#include "mesh.h"
#include "enumerations.h"
#include "constants.h"
// ***********************************************************************
// P2 element, conforming, 2D
// ***********************************************************************

static INT C_T_P2_2D_dof[3] = {1,1,0};
static DOUBLE C_T_P2_2D_nodal_points[12] = {0.0,0.0, 1.0,0.0, 0.0,1.0, 0.5,0.5, 0.0,0.5, 0.5,0.0 }; 
static INT C_T_P2_2D_Num_Bas = 6;
static INT C_T_P2_2D_Value_Dim = 1;
static INT C_T_P2_2D_Inter_Dim = 1;
static INT C_T_P2_2D_Polydeg = 2;
static bool C_T_P2_2D_IsOnlyDependOnRefCoord = 1;
static INT C_T_P2_2D_Accuracy = 2;
static MAPTYPE C_T_P2_2D_Maptype = Affine;

static void C_T_P2_2D_InterFun(DOUBLE RefCoord[], INT dim, DOUBLE *values)
{ 
  values[0] = 1-3*RefCoord[1]+2*RefCoord[1]*RefCoord[1]-3*RefCoord[0]+4*RefCoord[0]*RefCoord[1]+2*RefCoord[0]*RefCoord[0]; 
  values[1] = -RefCoord[0]+2*RefCoord[0]*RefCoord[0];
  values[2] = -RefCoord[1]+2*RefCoord[1]*RefCoord[1];  
  values[3] = 4*RefCoord[0]*RefCoord[1];
  values[4] = 4*RefCoord[1]-4*RefCoord[1]*RefCoord[1]-4*RefCoord[0]*RefCoord[1];
  values[5] = 4*RefCoord[0]-4*RefCoord[0]*RefCoord[1]-4*RefCoord[0]*RefCoord[0];
}

// base function values
static void C_T_P2_2D_D00(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ 
  //v =[1-3*x(:,2)+2*x(:,2).^2-3*x(:,1)+4*x(:,1).*x(:,2)+2*x(:,1).^2, 
  // -x(:,1)+2*x(:,1).^2,  -x(:,2)+2*x(:,2).^2,...
  //    4*x(:,1).*x(:,2), 4*x(:,2)-4*x(:,2).^2-4*x(:,1).*x(:,2),  4*x(:,1)-4*x(:,1).*x(:,2)-4*x(:,1).^2];
  values[0] = 1-3*RefCoord[1]+2*RefCoord[1]*RefCoord[1]-3*RefCoord[0]+4*RefCoord[0]*RefCoord[1]+2*RefCoord[0]*RefCoord[0]; 
  values[1] = -RefCoord[0]+2*RefCoord[0]*RefCoord[0];
  values[2] = -RefCoord[1]+2*RefCoord[1]*RefCoord[1];  
  values[3] = 4*RefCoord[0]*RefCoord[1];
  values[4] = 4*RefCoord[1]-4*RefCoord[1]*RefCoord[1]-4*RefCoord[0]*RefCoord[1];
  values[5] = 4*RefCoord[0]-4*RefCoord[0]*RefCoord[1]-4*RefCoord[0]*RefCoord[0];
}

// values of the derivatives in RefCoord[0] direction
static void C_T_P2_2D_D10(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ 
  //vx =[ -3+4*x(:,2)+4*x(:,1),  -1+4*x(:,1),  zeros(dim,1),    4*x(:,2),       -4*x(:,2),    4-4*x(:,2)-8*x(:,1)];
  values[0] = -3.0+4.0*RefCoord[1]+4.0*RefCoord[0];
  values[1] = -1.0+4.0*RefCoord[0];
  values[2] = 0.0;
  values[3] =  4.0*RefCoord[1];
  values[4] = -4.0*RefCoord[1];
  values[5] =  4.0-4.0*RefCoord[1]-8.0*RefCoord[0];  
}

// values of the derivatives in RefCoord[1] direction
static void C_T_P2_2D_D01(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{ 
  //-3+4*x(:,2)+4*x(:,1),   zeros(dim,1), -1+4*x(:,2),    4*x(:,1),  4-8*x(:,2)-4*x(:,1), -4*x(:,1)
  values[0] = -3+4*RefCoord[1]+4*RefCoord[0];
  values[1] = 0;
  values[2] = -1+4*RefCoord[1];
  values[3] = 4*RefCoord[0];
  values[4] = 4-8*RefCoord[1]-4*RefCoord[0];
  values[5] = -4*RefCoord[0];
}
// values of the derivatives in RefCoord[0]-RefCoord[0]  direction
static void C_T_P2_2D_D20(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
  values[0] = 4;
  values[1] = 4;
  values[2] = 0;
  values[3] = 0;
  values[4] = 0;
  values[5] = -8;
}
// values of the derivatives in RefCoord[0]-RefCoord[1] direction
static void C_T_P2_2D_D11(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
  values[0] = 4;
  values[1] = 0;
  values[2] = 0;
  values[3] = 4;
  values[4] = -4;
  values[5] = -4;
}
// values of the derivatives in RefCoord[1]-RefCoord[1] direction
static void C_T_P2_2D_D02(ELEMENT *Elem, DOUBLE Coord[], DOUBLE RefCoord[], DOUBLE *values)
{
  values[0] =  4.0;
  values[1] =  0.0;
  values[2] =  4.0;
  values[3] =  0.0;
  values[4] = -8.0;
  values[5] =  0.0;  
}
// values of the derivatives in RefCoord[1]-RefCoord[1] direction
static void C_T_P2_2D_Nodal(ELEMENT * elem, FUNCTIONVEC *fun, INT dim, DOUBLE* values)
{  
  //INT dim = 1;
  DOUBLE coord[2];
  INT i;
  INT j;
  INT k;
  for(i=0;i<3;i++)
  {
    // dof on the verts 
    coord[0] = elem->Vert_X[i]; 
    coord[1] = elem->Vert_Y[i];
    fun(coord, dim, values+i*dim);
    //dof on the midpostatic ints of the edges
    j = (i + 1) % 3;
    k = (i + 2) % 3;
    coord[0] = 0.5*(elem->Vert_X[j]+elem->Vert_X[k]); 
    coord[1] = 0.5*(elem->Vert_Y[j]+elem->Vert_Y[k]);
    fun(coord, dim, values+(i+3)*dim);
  }
  //printf("Use the P2 nodalfunction, dim=%d!\n",dim);
 }
#endif